MRI offers wonderful soft-tissue resolution and various tissue contrast depending on the sequence involved. It is necessary to understand a bit of its principle before we talk about intensity normalization.
Weighted-sequence
T1-weighted (T1w) and T2-weighted (T2w) are the most commonly acquired sequence using MRI. The voxel intensity of T1- and T2w images typically represents the time required by the spin to resume to the main-field alignment and to archive transverse plan relaxation, respectively.
Quick intro - Bloch equations for nuclear magnetization
Let the field experienced be only the main field \(\textbf{B}(t)=B_0\hat{z}\) such that the cross product term becomes \(-M_x(t)B_0\hat{y} + M_y(t)B_0\hat{x}\) and the equation becomes:
Solving this numerically is inefficient but simple. As a quick example, consider an initial magnitization of \(\textbf{M}(0) = (1, 1, 1)\), with \(\gamma = 1, T_1 = 4, T_2 = 8\), subjected to a constant main B-field of \(\textbf{B}(t) = \hat{z}\):
import matplotlib.pyplot as pltimport numpy as npfrom numpy.linalg import normfrom IPython.display import*steps =1000# number of stepsT =100.# simulation lengtht = np.linspace(0, T, steps) # time axisdt =t[1] - t[0] # time resolutionB = np.array([0, 0, 1.]) # Main B-fieldM_0 = np.array([0., 0., 1.]) # equilibrium spingamma =1T1 =4# These are the weighting that dictates how the spin relax to z-axisT2 =8# These are the weighting that dictates how the spin relax to z-axis and transversely in xy plane# Starting simulationpath = []M = np.array([1., 1., 1.])M /= norm(M)for i inrange(steps): bx, by, bz = B # ignore bx and by mx, my, mz = M dMx = gamma * my * bz - mx / T2 dMy =-gamma * mx * bz - my / T2 dMz =- (mz - norm(M_0)) / T1 M = M + np.array([dMx, dMy, dMz]) * dt path.append(np.copy(M))path = np.stack(path)# Plottingfig, ax = plt.subplots(1, 1, figsize=(5, 3))ax.plot(t, path[:, 0], label="$M_x$")ax.plot(t, path[:, 1], label="$M_y$")ax.plot(t, path[:, 2], label="$M_z$")ax.set_xlabel("Time (s)")ax.set_ylabel("Magnetization")plt.legend()plt.show()
Now while you see one arrow in the animation, there are many protons within the scan FOV and its impossible to considered each of them individually. Therefore, the term